!Vérification des routines d'intégration
module fonctions
implicit none
integer, parameter :: points = 10000
contains

!------------------------------------------------------------------------
!Méthode de simpson pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine simptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
double precision :: wgt
somme = y(0)
somme = somme + y(n)
wgt = 2.d0
if (n==0) then
	somme = 0.d0
end if
if (n > 1) then
	do i=1,n-1
	  if (wgt == 2.d0) then
	    wgt = 4.d0
	  else
	    wgt = 2.d0
	  end if
	  somme = somme + wgt*y(i)
	end do
	somme = h*somme/3.d0
else
	somme = h*somme/2.d0
end if
return
end subroutine simptab

!------------------------------------------------------------------------
!Méthode du trapèze pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine traptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n > 1) then
	do i=1,n-1
	  somme = somme + 2.d0*y(i)
	end do
end if
somme = h*somme/2.d0
return
end subroutine traptab


end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program annsimtest
use fonctions
implicit none
double precision, dimension(0:points) :: tableau
double precision :: xi,xf,h,somme
integer :: i

xi=0.d0
xf=100.d0
h = (xf-xi)/dble(points)
do i=0,points
	tableau(i) = dexp(-h*dble(i))
end do
call simptab(tableau,h,points,somme)
write(*,*) somme
call traptab(tableau,h,points,somme)
write(*,*) somme
write(*,*) "On voit que la formule de Simpson (première) est bien meilleure que celle des trapèzes (deuxième)."
read(*,*)

end program annsimtest